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Using exact continuous quantum Monte Carlo techniques, we study the zero and finite temperature 
properties of a system of harmonically trapped one dimensional spin 1 /2 fermions with short range 
interactions. Motivated by experimental searches for modulated Fulde-Ferrel-Larkin-Ovchinikov 
states, we systematically examine the impact of a spin imbalance on the density profiles. We quantify 
the accuracy of the Thomas- Fermi approximation, finding that for sufficiently large particle numbers 
(N > 100) it quantitatively reproduces most features of the exact density profile. The Thomas-Fermi 
approximation fails to capture small Friedel-like spin and density oscillations and overestimates the 
size of the fully paired region in the outer shell of the trap. Based on our results, we suggest a 
range of experimentally tunable parameters to maximize the visibility of the double shell structure 
of the system and the Fulde-Ferrel-Larkin-Ovchinikov state in the one dimensional harmonic trap. 
Furthermore, we analyze the fingerprints of the attractive contact interactions in the features of the 
momentum and pair momentum distributions. 

PACS numbers: 03.75.Ss,71.10.Pm,02.70.Ss 
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I. INTRODUCTION 

Mapping the phase diagram of fermions with attractive 
interactions and spin imbalance is a challenging problem, 
both from the experimental and the theoretical side. For 
an unpolarizcd system, an effective attractive interac- 
tion leads to the formation of Cooper pairs carrying zero 
momentum^ However in the presence of spin imbalance, 
where the different spin components have their Fermi en- 
ergies shifted from one another, the pairing mechanism 
changes. Depending on parameters, different theories ap- 
ply. The Sarma or breach-pair theory, which is stable 
in the presence of long range interactions, suggests a co- 
herent superposition of Bardeen-Cooper-Schrieffer (BCS) 
pairs and normal fermions,— i 3 -^ 5 - while in other regimes 
one finds Fermi surface deformations^ 7 - or a coexistence 
of normal fluid and superfluid^ One of the most inter- 
esting possibilities is the Fulde-Ferrel-Larkin-Ovchinikov 
(FFLO) states , 9 ' 10 ! 11 which reconcile spin imbalance and 
superconductivity by pairing particles whose momentum 
does not sum to zero. For example, in the FF state, f- 
spin fermions with momentum k are paired with |-spin 
fermions with momentum —k + q, resulting in a Cooper 
pair with momentum q. In the LO theory, the Cooper 
pairs carry no net current, but are in superpositions of 
momentum states, leading to an order parameter that 
oscillates in space; the excess spin polarization resides in 
the nodes of the order parameter. Analogs of the lat- 
ter state are found in one dimension (ID). We will use 
an exact quantum Monte-Carlo procedure to study a gas 
of spin imbalanced harmonically trapped ID fermions, 
with a goal of identifying the signatures of FFLO physics 
which will be seen in cold atom experiments. 

Although not completely transparent, the best evi- 



dence for the FFLO state in solid-state systems comes 
from layered organic superconductors^ and heavy- 
fermion materials . 13 ' 14 However, recent progress in op- 
tical lattice experiments with cold atoms has opened the 
avenue to study spin imbalanced fermions with attrac- 
tive interactions in a clean and direct wa y 15 ' 16 ' 17 ' 18 ' 19 ' 20 
Unfortunately, for a bulk three dimensional gas of atoms 
with short range interactions, the FFLO state is found 
for a very narrow range of chemical potentials^ imply- 
ing that only a small volume of a trapped atomic cloud 
would be in this state. Hence attention has been shifting 
to one dimension, where an analog of the FFLO state has 
a very large region of stability^ 

Since there is no true off-diagonal long-range order in 
a strict ID system, the ID FFLO state is characterized 
by an algebraic quasi-lD order. Although the thermo- 
dynamics of a uniform ID Fermi gas with contact in- 
teractions can be found exactly via the Bethe ansatz 
at zero ' and finite temperature , 25 ' 26 the characteriza- 
tion of the ground state with spin imbalance is not easy: 
many correlation functions are unknown due to the com- 
plicated structure of the Bethe ansatz. Other properties 
are well known: for example, Yang firmly established the 
existence of the FFLO state through a bosonization tech- 
nique, showing that it emerges via a continuous transition 
from the BCS stated Further, numerical methods have 
shown that the FFLO is the ground state of the Hamilto- 
nian on a lattice as soon as the system is polarize d 28 ' 29 ' 30 . 

Experimentally, an array of quasi one dimensional 
tubes can be created at the intersection of two orthogonal 
standing waves generated by laser beamsi 31 ' 32 ' 33 Along 
each tube there would be no lattice, but the particles 
will be confined by a longitudinal harmonic potential, 
coming from the spatial profile of the lattice lasers, and 
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possibly from some additional potential. Recently Orsc 34 
and Liu at al^ studied the properties of an atomic cloud 
in such a ID harmonic trap. By using a Thomas Fermi 
(TF) approximation based on the Bethe ansatz solution 
of the homogeneous Hamiltonian, they found that in the 
trap the system phase separates. The inner shell of the 
tube is always a partially polarized FFLO state, while 
the outer shell can be either a fully paired BCS or a fully 
polarized normal state. Our calculations are designed to 
go beyond the approximations inherent in those calcula- 
tions, and accurately quantify the properties of these ID 
clouds. 

Other researchers have also attempted to investi- 
gate the accuracy of the TF approximation. For ex- 
ample, Xianlong and Asgari developed an exchange- 
correlation functional to perform density functional the- 
ory (DFT) calculations with the local density approxi- 
mation (LDA). 36 While revealing, those calculations are 
only approximate. In a different direction, several re- 
searchers have used the density matrix renormalization 
group (DRMG) to study a related system, where there 
is a ID lattice along the tubej 37 ' 38 i 39 One needs to be 
cautious about applying those results to the continu- 
ous system as the underlying lattice in the model could 
introduce effects of incommensurability in the density 
profilesj 40 ' 41 The continuous limit is difficult to access 
due to the large number of sites required in the DMRG 
evaluation for an extremely dilute system. 

We will also consider how temperature affects the den- 
sity profiles. Since the superfluid critical temperature is 
zero in ID, thermal effects could be quite dramatic, and 
may not be captured by mean field theories^ Experi- 
ments will inevitably be performed at finite temperature, 
and it is imperative to know to what extent finite temper- 
ature will obscure the clear delineation of phases which 
are expected in a trapped zero temperature gas. 

In this paper, we present exact results for ID fermions 
subject to a zero-range attractive interaction and a lon- 
gitudinal harmonic confinement. We resolve the de- 
pendence on the spin imbalance, effective coupling, and 
temperature, within a continuous quantum Monte Carlo 
(QMC) framework. We use diffusion Monte Carlo 
(DMC)^ 3 - for zero temperature ground state calculations 
and path integral Monte Carlo (PIMC)^ for finite tem- 
perature, which are ideal numerical tools to find the prop- 
erties of ID Hamiltonians, as they do not suffer from the 
"sign problem" which affects QMC simulations in higher 
dimensions^ In ID they give unbiased energies, with 
their accuracy only limited by statistical noise. We cor- 
rect other potential errors in the DMC algorithm, such 
as the time step r, population bias, and mixed estimate 
errors, by appropriate extrapolation and forward walk- 
ing techniques^ The only systematic bias in the PIMC 
algorithm is the time step, which we control by working 
with sufficiently small r. 

The paper is organized as follows. In Sec. [IT] we outline 
the model Hamiltonian used in our calculations, and the 
experimental parameters of the system. Sec. Mil includes 



a description of the QMC methods with an emphasis on 
the features required to simulate ID fermions interact- 
ing with a zero-range potential. In Sec. IIVI we present 
our results and compare them with the TF findings. In 
particular, we focus our attention to the finite size and 
temperature effects on the phase boundaries. We con- 
clude in Sec. [V] with some highlights and remarks. 

II. MODEL 

The Hamiltonian which describes ID fermions in a trap 
and interacting with an attractive zero-range potential is 

h 2 N d 2 nT Nl 1 N 

i=X 1 i—1 j—1 i—1 

(!) 

where N — + is the total number of particles with 
mass to, <?(> 0) is the strength of the contact interaction, 
and Lo z is the harmonic frequency of the longitudinal trap. 
The effective ID Hamiltonian in Eq.Q]is a very good rep- 
resentation of a single tube in the optical array provided 
that Nhoj z <C huj±_ and fc^T <C tko±, where lo± is the 
frequency of the transverse harmonic confinement. The 
above requirements are met in the experimental setup 
by properly tuning the height of the lattice depth, which 
sets the value of co±. In such a ID geometry with just 
one open channel, the interactions can be modeled by as 
V(x) — ~gS(x) with g given b y 47 ' 48 

2H 2 a 3D 1 
9 = — -j / 7, ( 2 ) 

w here a 3 p is the 3D s-wave scattering length, a± = 
yJKjmojZ is the transverse oscillator length, and A — 
— £(1/2)/a/2 ~ 1.0326, where £ is the Riemann zeta func- 
tion. Here and throughout the paper we have used a sign 
convention where g > is attractive. An ongoing exper- 
iment at Rice University 4 ^ is using N ~ 100 — 300 of 6 Li 
atoms per tube, the actual value depending on the posi- 
tion of the tube in the array. Typical values for the opti- 
cal lattice arc u±_ = 2ir 156 kHz, uj z — 2tt 200 Hz, while 
g will be tuned by changing a 3 u- The typical tempera- 
tures are around T ~ 0.05 T 1 F , with T% = N a hu; z /k B the 
Fermi temperature of fermions with spin a. Throughout 
the paper we will use m = % = u> z = ks = 1 ; i-e. the 
length will be measured in units of the harmonic oscilla- 
tor length ^h/mujz, and the energy and temperature in 
units of the level spacing htu z . In these units the energy 
for g = is E = (N 2 + Nf)/2 and the binding energy of 
a pair is g 2 /4. 

III. METHODS 

The ground state QMC techniques used in this work 
are described in Subsec. IIII A[ while the finite tempera- 
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ture PIMC method is reported in Subsec. IIII Bl In or- 
der to compare our results with the local density ap- 
proximation, we computed also the TF density profile 
(n(x) = rv(x) + n^(x)) and spin density profile (s(x) — 
n^{x) — n^(x)) of the system, which are solutions of: 



=U, (3) 



where \i a is the chemical potential of the species with 
spin tr, and /Zhom.o- is the local chemical potential found 
from the Gaudin exact solutio n 23 i 24 of the homogeneous 
model. The Thomas-Fermi results have been widely dis- 
cussed elsewhere, and we refer the reader to previous 
works 22 ' 34 ! 35 ' 36 for complete details. Here, we mention 
that Eq. [3] can be rescaled into dimensionless variables 
with coupling strength controlled by 4N/g 2 . Therefore, 
the effective coupling constant of the trapped system 
depends also on the number of particles and scales as 
9/VN. 



A. Ground state QMC 

To study the ground state properties of the system, 
we employed variational and diffusion quantum Monte 
Carlo methods. Their basic ingredient is the trial wave 
function ^t, & good approximation to the ground state. 
We used the form: 



* T = £> T D l exp(-fZ), 



(4) 



where U is a bosonic function and D a is the antisym- 
metrized product of Hermite polynomials of i-th order 
Hi for particles with spin a. D" can be written as deter- 
minant of a Van der Monde matrix: 

D a (xi, . . . ,x N <?) = det(Hi(xj)) = {xi - Xj). 

l<i<j<N" 

(5) 

The product in the Eq.[5]is convenient, as it can be evalu- 
ated in order iV 2 operations. The Hermite basis set yields 
an exact variational wavefunction in the limit of vanish- 
ing inter-particle interactions. The Slater term 
determines the nodes of the variational wave function. 
In ID, the exact nodes of the ground state are defined by 
the coalescence conditions between two like-spin particles 
(i.e. Xi = x.j with di = <r,j). This result holds even for 
an interacting system, and consequently there is no sign 
problem in the diffusion Monte Carlo calculations^ 5 - 

The bosonic function U is real and symmetric under 
permutation of two fermions, and is the sum of one-body 
{Uib), two-body (U 2 b), and three-body {U 3 b) correlations. 
The one-body part reads: 



N 



(0) 



where Ui is the spin of the i-th particle, and f(x) is 
written in a Pade form which goes as x 1 jl at large x 



so to assure the right asymptotic behavior given by the 
harmonic confinement. In the absence of inter-particle 
interactions, f a {x) = x 2 /2 is the exact solution. 
The two-body (Jastrow) factor is defined as: 



U 2b = 2^ U°*°i ( Xi - Xj ) , (7) 

l<i<j<N 

with u aiGi pair functions chosen to fulfill the cusp condi- 
tions set by the contact interactions of the Hamiltonian 
in Eq. [1] In the case of two unlike-spin particles inter- 
acting with a delta function potential, the exact solution 
is given by ^{x} ,x^) = exp(— — % \)', note the cusp 
at x^ — x- 1 . The exact many-body function will have the 
same cusp between atoms with different spin. We use the 
following pair functions: 



,tl 



(x) = -— exp(-A|a;|) 



2A 



u aa {x) 



exp(-/3 aa \x\) 



(8) 

^{U},(9) 



' ff "l + cxp(-2/3 (T(T |s|)' 

where A, a aa , and [3 aa are variational parameters, and g 
is the strength of the contact interaction. 

The confining potential leads to inhomogeneities that 
make it convenient to include an inhomogeneous Jastrow 
function: 



^non-homo = £ ( Xi - Xj ) (x i+Xj ), (10) 



l<i<j<N 

where the pair function u (Ti(7j is modulated by the func- 
tion / <Ti<T J which depends on the center of mass of the 
particles Xi and Xj. The analytic form of u aa is the 
same as in Eq. [9] also for a ^ er', while f aa is chosen to 
be: 

exp(-S a(T 'x) 



J™ (x) = l + la 



(11) 



1 + eXp(-2<5 c r C r'2!) ' 

where 7 CT(T ' controls the inhomogeneity. For 70-0-' = one 
recovers the homogeneous case. 

Finally, we employed a three-body term: 

N 

U 3b = u^ixi-Xj) u^-ixj-Xk), (12) 

where also in this case the functional form of u is the 
same as in Eq. [SJ We found this term important in the 
intermediate and strong coupling regime, when the local 
energy is dominated by the inter-particle potential. 

Our variational ansatz includes about 30 parameters 
{pi}, optimized by means of the stochastic reconfigura- 
tion algorithm ; 50 ' 51 which minimizes the total energy Et 
of the wave function with an iterative procedure based on 
the value of the forces, —dEx/dpi, acting on the param- 
eters at each step. The optimal variational wave function 
is then projected to the ground state by using the diffu- 
sion Monte Carlo method (DMC)^ 3 -, which provides an 
unbiased estimate of the ground state energy. The phys- 
ical observables, such as the charge and spin density pro- 
files, are computed with the forward walking propagation 
of the mixed DMC distribution^ 6 .^ 
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B. Finite temperature QMC 

The finite temperature calculations are performed by 
using the path integral Monte Carlo (PIMC) method^ 
which is based on sampling the thermal density ma- 
trix p{R,R'\[3), with (3 = 1/fcgT, T the temperature, 
by means of a mapping of the quantum system into 
a classical system of polymers or Feynman 53 "paths" . 
R = {x\, . . . , xn} and R' = {yi, . . . , ijn} are all-particle 
configurations. During the Monte Carlo random walk the 
polymers are distributed according to the probability 

exp ( - J2 Sm) , (13) 

\ m=l / 

where the integration from to j3 has been discretized 
into M "time slices" of length t = /3/M, S m = 
— log [p(i? m _i, Rm\ t)] is the action, and Rq = Rm- An 
accurate approximation of p(R, R '; r) for the system un- 
der consideration is important to make the PIMC calcula- 
tions feasible and efficient. In the case of the Hamiltonian 
in Eq. [TJ p(R,R';r) can be written using a Trotter ex- 
pansion as a product of four factors: the non-interacting 
(po), the harmonic (pharm), the interacting (pint), and the 
fcrmionic (pFermi)- The non-interacting part corresponds 
to the propagator for free particles 

Po (R, R; t) = (2 7 rr)- JV / 2 exp {-(R - R') 2 /2r} . (14) 

We use the primitive approximation^ to include the ef- 
fect of the harmonic external potential: 

Pharm(i?, R' '] t) = exp {-1 (R 2 + R' 2 ) } . (15) 

To account for the interactions we make use of the exact 
action for two particles in free space interacting with an 
attractive contact potential^. In terms of the relative 
distances of the pair (x rc \ = x% — x-2, y Ic i = Vi — y-i), the 
resulting propagator is: 

ps(x re l, J/rei; t) = 1 + y/vS exp {s(s - 2r) 

+z 2 (x lcl - y Icl ) 2 } erfc (r - s) (16) 

where s = \yjr, z = 1/V4r, r = z (\x Te \\ + \y IC \\). From 
this propagator, one obtains the many-body action by 
means of the pair-product approximation^ 

Pint(i?,i?';r) = JJJJp5(xi -Xj,yi-yf,T), (17) 
i=ii=i 

where the product runs over all unlike-spin pairs. This, 
as the other approximations, becomes exact as r — > 0. 

To include Fermi statistics into the PIMC method^ 
the paths are constrained by the nodal regions, which 
act like barriers. Since the nodes are exactly known in 
ID, and they do not depend on /3, the restricted path in- 
tegral formalism gives the exact result. An approximate 



fermion action is given by the image approximation: 

p Fcrmi (i?, R ; r) = I M 1 - exp i — | 

(18) 

for each spin sector. It includes the nodes between neigh- 
boring particles as zero's of the density matrix, and be- 
comes exact in the r = limit. We assume that the like 
spins are ordered Xi-\ < Xi for 1 < i < N a . 

The probability distribution in Eq.[l3l based on the the 
density matrix pi n t(-R, R'\ t) defined above, is sampled via 
a generalized Metropolis algorithm. The paths are dis- 
placed with two types of moves: a multi- level bisection 
(level two and three have been used here) and a global 
move, which attempts to displace an entire polymer at 
once. When the convergence is reached, the thermal ex- 
pectation values are computed as described in Ref. |4J. 

IV. RESULTS 

The QMC results shown in this section characterize 
the ground state and finite temperature properties of the 
system. We compute the density n(x), spin density s(x), 
and local polarization (p(x) = s(x)/n(x)) profiles. The 
latter quantity is extremely useful, as it locates the paired 
(p = 0), partially polarized (0 < p < 1), and fully po- 
larized (p = 1) shells. For sufficiently high densities, 
the center of the cloud is always partially polarized, but 
the edge may be paired or fully polarized. The Thomas- 
Fermi approximation yields the phase diagram shown in 
Fig.[TJ which shows the expected boundary between these 
two regimes. We decided to predominantly study the sys- 
tem with N = 200, since it is roughly the average number 
of fermions per tube in typical experiments. As marked 
on the phase diagram, we present results for a range of 
interaction strengths and polarizations. In our numerical 
profiles we are particularly interested in how the visibility 
of the interface between different regions is influenced by 
finite particle number and finite temperature. By ana- 
lyzing the experimental data within a Thomas- Fermi ap- 
proximations one could use the locations of these bound- 
aries to map out the phase diagram. We find that for 
N > 100 and T < 0.05T F this Thomas-Fermi analysis is 
reasonable, but for fewer particles or higher temperatures 
one would need to use a much more sophisticated analy- 
sis of the experimental data to learn about the bulk zero 
temperature phase diagram. We report our results for 
the ground state and the finite temperature calculations 
in Sees. IIV Al and HV Bl respectively. 

A. Ground state 

In Fig. [2] we report on the density, spin den- 
sity, and local polarizations of N = 200 particles 
with various polarizations P and interaction strengths 
g. From the TF results in Fig. [TJ one expects 
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FIG. 1: (Color online) Thomas-Fermi phase diagram in 
the polarization-coupling plane for ID harmonically trapped 
cloud of two-component fermions, interacting through a con- 
tact potential V(xi — Xj) = —g8(xi — Xj). The number of 
particles is N, the polarization P, and we use units where the 
harmonic confinement energy and length scales are set equal 
to unity. The gray thick line separates the phase region where 
the outer shell of the cloud is fully paired (S) from the region 
with fully polarized outer shell (P). The points locate our 
QMC simulations on the TF phase diagram. Triangles (red 
outline), squares (blue outline), and diamonds (green outline) 
refer to DMC data reported in Fig. [2] panels (a) , (b) , and (c) 
respectively. The two circles at P = 0.10 depict systems with 
the same g = 20 but different sizes (N = 60 and 20), displayed 
in Fig. [3] The circle at P = 0.04 represents the system with 
N = 200 and g — 16, whose ground-state DMC polarization 
is reported in Fig. [4] At the same polarization, the trian- 
gle, circle and square points have been studied also by PIMC, 
which yielded temperature dependent polarizations reported 
in Figs. [SHU and [10] respectively. 



that the edge of the cloud will be fully polarized 
for (g,P) = (8, 0.2), (8, 0.15), (8, 0.1), (20, 0.2), (20, 0.15). 
Similarly, the edge should be fully paired for (g, P) = 
(20, 0.04), (50, 0.1), (50, 0.04). The others: (g,P) = 
(8, 0.04), (20, 0.1), (50, 0.15) are near-critical, and one ex- 
pects that the FFLO state should extend nearly to the 
edge of the cloud. The results in Fig. [2] are consistent 
with these expectations. For the TF near-critical points, 
we found that all of them feature fully polarized wings. 
When the outer shell of the cloud is polarized, the spin 
density has a peak at the edge, which is replaced by a 
rapid fall-off when the edge is fully paired. Detecting the 
fully paired shell will be easiest at strong coupling, where 
the size of the outer paired region can be large. 

Notice that the g = 50 spin polarization plotted in 
Fig.[2jc) is much more noisy than for weaker interactions. 
Indeed, from the computational point of view a problem 
of ergodicity starts affecting the Monte Carlo (MC) sam- 
pling. As the coupling increases, one needs to run longer 
to equilibrate the sample and have good statistics in the 
expectation values, because sometimes the evolution gets 
stuck in a "persistent" configuration. In this regime, it is 
difficult for unpaired fermions to cross a strongly bound 



pair due to the Pauli principle. We note that it is pos- 
sible that the same physics will affect also the dynamics 
of the 6 Li atoms in the experiment. The physical origin 
of the ergodicity problem in the MC evolution, namely 
the formation of strongly bound pairs, could also cause a 
slow equilibration in the experimental setup. A coupling 
between g = 20 and g = 50 is a good balance between 
the visibility of the fully paired region and the ergodicity 
of the pairs in the tubes. 

From the analysis presented above, it seems that the 
TF mean field phase diagram is a good approximation 
for the couplings and polarizations studied so far. It is 
clear that the TF approximation is not able to reproduce 
the short length-scale charge and spin fluctuations, since 
the local density approximation assumes a slow vary- 
ing density profiles, as already pointed out in previous 
studies . 36 ! 56 However, the TF boundaries of the phase 
separated states are in a quite good agreement with our 
exact QMC calculations. In order to better understand 
this agreement, we study the dependence of the local po- 
larization on the number of particles N and coupling g, 
and compare our QMC results with the TF predictions. 
The local polarization for different system sizes (N = 20, 
60, 200) with g = 20 and P = 10% is plotted in Fig. [3] 
At N — 20 and 60, there are qualitative disagreements 
between the QMC and TF methods, which disappears 
at N = 200. Indeed for small systems our QMC cal- 
culations give partially polarized wings, while they are 
fully paired in the TF approximation. This can be due 
to a proximity effect not accounted for in the completely 
local TF approximation. We speculate that the FFLO 
state partially extends into the wings, making them par- 
tially polarized. We cannot confirm this picture since we 
did not study how the FFLO order parameter depends 
on position. For N — 200 the wings become fully po- 
larized, as the effective coupling is reduced, and the po- 
larization profile turns out to be in agreement with the 
TF results, although the transition from the FFLO to 
the fully polarized state is much smoother in our QMC 
calculations. The boundaries between phase separated 
states are smeared out over a length-scale of order the 
atomic spacing. For small systems the very concept of 
phase separation begins to break down. 

We check the agreement between the QMC and TF 
methods for large number of particles (N = 200) in a 
more systematic way by studying the dependence of the 
polarization on g with P — 4%. These results are plotted 
in Fig. g] for g = 8, 16, 20, and 50. The TF polarization 
profile follows quite closely the DMC points. Both the- 
ories predict a "transition" from a fully polarized shell 
to a fully paired one between g — 8 and 16. Apart from 
fluctuations and statistical noise on the points, the only 
significant difference in the comparison is the size of the 
fully paired region which is smaller in the QMC results, 
where the FFLO state seems to extend farther from the 
center of the trap. This is a common feature of our calcu- 
lations. The TF approximation always overestimates the 
size of the fully paired shell. Apparently, this is a draw- 
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FIG. 2: (Color online) Density n, spin density s and local 
polarization p = s/n profiles for a cloud with N = 200. The 
panel (a) shows the QMC results for g = 8 and four different 
polarizations P, while panels (b) and (c) display the QMC 
density profiles for g = 20 and 50 respectively, with the same 
set of polarizations. The spin density is reported at the right 
side of the plot, while the density is on the left. The inset 
shows the polarization at the edge of the cloud, where the 
crossover from the FFLO core to the outer shell occurs. At 
g = 8 the edge of the cloud is a completely polarized normal 
gas for all polarizations, whereas for g = 20 and the lowest 
polarization P the edge is in the fully paired state. At even 
stronger coupling (g = 50), both P = 4% and P = 10% give 
fully paired outer shells. Finally, note how the strength of 
the interaction shrinks the cloud, while up to P — 20% the 
polarization has a little impact on its spread. At g — 50 the 
equilibration of the system begins to break down in the MC 
sampling, as particles are unable to easily cross each other, 
being the pairs strongly bound. This breakdown accentuates 
the MC noise in the results. 



back of the local density approximation, which makes the 
effective attractive interaction stronger. 57 




4 6 
distance from the center 

FIG. 3: Local polarization profiles for g — 20, P = 10%, 
and different number of particles N. The TF polarization is 
drawn (thick black line) for comparison. The TF approxima- 
tion gives fully paired wings with zero polarization for N = 20 
and 60, while the QMC polarization is always positive. The 
agreement between QMC and TF is better for the largest sys- 
tem (N = 200), where the outer shell is fully polarized in both 
methods. 
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FIG. 4: Local polarization profiles for N = 200, P = 4%, 
at various couplings g. Also the TF polarization is drawn 
(thick black line) for comparison. The QMC results show a 
crossover to fully paired wings for g > 16, which is much 
smoother than the corresponding TF transition. Moreover, 
the QMC simulations give a larger partially polarized shell in 
the center of the trap. 

We complete the comparison between the TF and 
QMC results by also computing the energetics of g = 20 
with different polarizations (P — 0%, 10%, 20%) and sys- 
tem sizes (N = 20, 200), reported in Tab. [fl The agree- 
ment in the ground state energies per particle between 
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the two methods is remarkable, even at strong coupling 
(small N). This means that the total energy is insen- 
sitive to the local differences seen in the shell structure 
of the cloud, interpreted as finite size or edge effects. It 
also highlights the importance of an accurate treatment 
of local correlation effects, which can be missed in the TF 
approximation, and lead to subtle changes in the phase 
boundaries, whereas the overall trend is correctly cap- 
tured by the mean field method, even at a quantitative 
level. 

TABLE I: Ground state energies per particle computed by 
means of DMC and TF, for g = 20 at different polarizations 
(P = 0%, 10%, and 20%) and system sizes (N = 20, and 
200). 



N 


P 


TF energy 


DMC energy 


20 


0% 


-47.265 


-47.270(1) 


20 


10% 


-42.550 


-42.538(3) 


20 


20% 


-37.587 


-37.585(2) 


200 


0% 


-35.731 


-35.478(10) 


200 


10% 


-13.991 


-13.983(3) 


200 


20% 


-7.914 


-7.899(3) 



In previous workS ) 28 i 29 ' i 38 it has pointed out that the 
pair momentum distribution function reveals fingerprints 
of the FFLO state in the center of the trap. The momen- 
tum and pair momentum distributions can be extracted 
in experiments from a statistical analysis of time-of-flight 
images of the expanding cloud of atoms when the confin- 
ing potential is removed. Here we study the momentum 
distribution function defined as 

n(k) =J dxdyp(x,y)e lHx - y \ (19) 

where p(x, y) is the one-body density matrix 

p(x, y) = J dx 2 ■ ■ ■ dx N ^(x, x 2 , • • • , x N )^(y, x 2 , ■ ■ ■ , x N ), 

and \& is the ground state normalized such that 

J dxp{x,x) = N. (21) 

In the actual DMC calculations, the density matrix in 
Eq. [201 is computed as follows 

p(x,y)*J2Ux-Xi) * T (''' ,y, ''\ ) , (22) 

where (• • • )ma indicates the MC averages over 
{xi, . . . , Xi, . . . , xn} sampled from the DMC mixed dis- 
tribution. Although the forward walking technique for 
non-local operators is possible, it has rarely been applied 
to the evaluation of the momentum distributions^ In 
our case, we checked that the density matrix p(x, y) com- 
puted at the variational MC level is close to the DMC 



mixed average (MA), due to the accuracy of the trial 
wave function ^>t- Therefore, our best estimate of p(x, y) 
is a mixed average with a small bias. In Fig. [5] we plot 
n(k) for different couplings and N = 200. We found that 
the long-range tail of the momentum distribution decays 
as a power law as soon as the attractive contact interac- 
tion is switched on, and it is fitted well by 1/fc 4 . This is 
the same behavior found by Minguzzi at al^ for the tail 
of the momentum distribution in a one dimensional gas 
of hard point-like bosons (Tonks gas) inside a harmonic 
trap. The power law decay is due to the cusp conditions 
of the wave function, a consequence of the attractive (or 
repulsive) contact interactions. The same decay has been 
shown recently by Santachiara and Calabrese for a one 
dimensional gas of impenetrable anyons^ As one can 
see in Fig. (SJ the momentum distribution does not dis- 
criminate between different phase separated states, but 
from the height of its tail one can get information on the 
effective interaction strength in the trap. 
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FIG. 5: (Color online) Momentum distribution for N = 200, 
P — 4%, at various couplings g, plotted in a log- log scale. 
Also the momentum distribution of the non interacting sys- 
tem is reported. The tails of the interacting n(k) follow a 
straight lines with a slope compatible with a 1/fc 4 decay. Mo- 
menta are measured in units of y/hmu) z , where lj z is the har- 
monic oscillator frequency along the tube. For comparison 
the characteristic momentum corresponding to the density at 
the center of the trap is fc^ = n/(2n) = 16.5 for g = 20, while 
fci? = 14.1 for the non- interacting system. 

On the other hand, the pair momentum distribution 
carries fingerprints of the FFLO state. We define the pair 
two-body density matrix as 



P P mr{x,y) = 




■ ■ ■ dxN^>(x' , x^, Xs, • • ■ ,xjy) x 



*(l/W,s 3 ,-" ,zjv), (23) 

where the pair {x^ , x^ } is chosen by selecting x^ and find- 
ing x^ as the position of the closest unlike-spin particle. 
x (y) is the center of mass of the {x' , x+} ({y, y^}) pair. 
This is an unambiguous and generic way to define the 
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pairs in a continuous system without introducing explic- 
itly a characteristic length. Once p pa h(x, y) is computed, 
the pair momentum distribution function n pa i r (fc) is eval- 
uated via Eq. [THl with p(x,y) replaced by p p& i r (x, y). 
Also the normalization is based on the same formula as 
in Eq. [211 but with N replaced by iVp a i r , he. the num- 
ber of pairs in the system. The results for n pa i r (£:) are 
presented in Fig. [5] for g — 20, N = 200, and various 
polarizations. A clearcut signature of the FFLO state 
is the presence of two peaks at k = ±\k F — kp\, which 
signals a pairing with non zero total momentum given 
by the difference between the Fermi surfaces of up and 
down spin particles. Since the trap breaks the transla- 
tional invariance, the Fermi momenta are approximated 
by k F — f^-j where L is the effective spread of the 
cloud. The location of the peaks in Fig. [6] follows closely 
the relation ±|fct — kL\. 




FIG. 6: (Color online) Pair momentum distribution (Fourier 
transform of the pair two-body density matrix in Eq. I23|l for 
g = 20, N = 200, at various polarizations P. Momenta are 
measured in units of \/hmuj z , where ui z is the harmonic os- 
cillator frequency along the tube. The peaks of the pair mo- 
mentum distribution are approximately given by i\k F — k F \, 
and are the signature of the FFLO state. 



Fig. [7] displays the n pa i r (fc) dependence on g at fixed 
polarization P = 4%. At weak coupling (<? — 8) the 
peaks are short, and in the range [— k F , kp] the pair mo- 
mentum distribution is almost flat. On the other hand, 
at g — 20 the two peaks at ±\k F — kp\ are taller, and 
the atoms in the cloud have a quite strong FFLO char- 
acter. In the intermediate case analyzed here, namely 
g = 16, the n pa i v (k) shows a fluctuating behavior, with 

peaks at ±\k F — k F \ but also at k = 0. By computing 
the Fourier transform (Eq. I19[) of p pa ,i v (x,y) constrained 
in the region with x,y € [— i? C utoff, -Rcutoff] (and plotted 
in the inset of Fig. [7]), we proved that the FFLO peaks 
come from the inner shell of the cloud, while the central 
peaks come from the wings. Indeed, around g = 16 the 
system undergoes the transition from fully polarized to 
fully paired (BCS) wings (see Fig. g]). Both the FFLO 



pairing in the center and the BCS pairing in the outer 
shell are weak. Going from the center to the edge, the 
crossover between the two gives rise to this fluctuating 
pattern in the pair momentum distribution. 
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FIG. 7: (Color online) Pair momentum distribution for P = 
4%, N = 200, at various couplings g. The functions have 
been rescaled such that the height of their peaks is around 1 
to make the comparison easier. The inset displays the pair 
momentum distribution for g — 16 computed from p pa i I (x,y) 

with X,y e [— iZcutaff) -Rcutoff]. 



B. Finite temperature 

We performed finite temperature PIMC calculations 
for g = 8, 16, and 20, with P = 4%, and N = 200. From 
the DMC analysis and the TF results we know that for 
g = 8 the phase separation is between an FFLO center 
and a fully polarized wing, while for g = 16 and 20 the 
outer shell is fully paired. Here we study how the phase 
separation evolves with the temperature. 

In Fig. [5] we plot the local polarization profile for 
g = 8. At a temperature above 0.065 T F , the fully po- 
larized shell becomes partially polarized with polariza- 
tion progressively reduced as the temperature is raised. 
At T ~ 0.40 T F the polarization is quite homogeneous 
throughout the cloud. 

The finite temperature results at g — 16 (20), reported 
in Fig. (jTDJ) , show a more complicated pattern. Up to a 
temperature of T ~ 0.06 T ] F (0.10 T F ), the polarization 
profile has a dip around x — 11.5 (10.5), which evolves 
into a fully paired region (with zero polarization) below 
T ~ 0.025 T F (0.035 rj,). When the dip is present, 
the polarization in the outer part of the cloud increases 
until a fully polarized state occurs at the very edge of 
the cloud. This esternal region involves only 0.1% of the 
atoms in the cloud, and so it will be hard to detect in 
the experiment. At higher temperatures, when the dip is 
gone, the fully polarized edge becomes unstable. For T ~ 
0.40 T F we have a situation analogous to the case with 
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g = 8, with a partially polarized state homogeneously 
spread over the cloud. 

A feature common of all couplings is a crossover tem- 
perature T c above which the phase separation disappears 
and the whole cloud becomes partially polarized. This 
temperature depends on the polarization, since it scales 
roughly as fv — ^ , whereas its dependence on the cou- 
pling is quite weak, at least until g — 20. We found 
T c ~ 0.065 Tp for g = 8, and T c ~ 0.10 T p for g = 20 
with P = 4%. In case of fully polarized wings, one can 
define another critical temperature based on the stability 
of the pairs at the edge of the cloud. The fermions disso- 
ciate into a partially polarized fluid if the temperature is 
raised above a threshold proportional to g 2 , as it depends 
on the binding energy of each pair. 

As in the DMC calculations, the PIMC simulations suf- 
fer from the ergodicity problem mentioned before. This 
is an issue of the strictly one dimensional Monte Carlo 
sampling, where the exchange between particles or pairs 
is suppressed at strong coupling. For this particular sys- 
tem, we found that this issue is less severe in the DMC 
simulations than in the PIMC ones. Improved moves or 
working in the grand canonical ensemble could alleviate 
or solve the problem. Here, however, we could afford sim- 
ulations with N = 200 at quite strong coupling, namely 
up to g = 50 in the DMC, and g = 20 in the PIMC 
framework. Therefore our analysis has not been limited, 
since we were able to study the intermediate-strong cou- 
pling regime, where the fully paired region shows up at 
not-so-small polarizations. 




distance from the center 

FIG. 9: (Color online) Local polarization profile for N = 
200, P — 4%, g — 16, at various temperatures T measured 
in units of Tiui z - For comparison, the central density in the 
trap corresponds to an energy ef = Ti 2 /2mk F ~ 259. The 
inset magnifies the region of the cloud around the polarization 
dip, where the pairing is the strongest. For temperatures 
above 2.5 the fully paired region starts breaking. When a 
pair is broken, the unpaired fermions feel the difference in the 
chemical potential between the up and down species, resulting 
in the majority of up fermions at the edge of the trap, and 
a fully polarized outer shell. For even higher temperatures, 
the cloud becomes partially polarized everywhere. The filled 
triangles indicate the polarization at the cloud radius which 
contains 99.9% of atoms. 
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FIG. 8: (Color online) Local polarization profile for N = 
200, P = 4%, g = 8, at various temperatures T measured 
in units of Tiuj z . For comparison, the central density in the 
trap corresponds to an energy ef = Ti 2 /2mk F ~ 229. At 
this coupling the outer shell is fully polarized only at the 
two lowest temperatures considered here, while it becomes 
partially polarized at a temperature in between T = 7 and 
T = 8. The filled triangles indicate the polarization at the 
cloud radius which contains 99.9% of atoms. The remaining 
0.1% (0.2 atoms) at the edge will be hard to detect in the 
experiment. 



V. CONCLUSIONS 

In this paper we have studied the properties of a system 
of ID trapped fermions interacting via an attractive con- 
tact potential, by performing unbiased DMC and PIMC 
simulations. We showed that the trap induces the system 
to phase separate, in accordance to previous calculations 
done with the TF approximation] 34 ! 35 We compared the 
TF solution with our exact results, by studying their de- 
pendence on the coupling and system size. We found 
that the DMC phase boundaries are in a good agreement 
with the TF approximation for large number of particles 
(N > 100) even at strong coupling. This is reasonable 
since the TF approach is applicable when N ^> 1, and 
the size of the cloud is large compared to the axial oscil- 
lator length. However, the local density approximation 
requires also a slow variation of the mean interparticle 
spacing, namely d x n(x) / n(x) 2 <C 1, condition which is 
not fulfilled by the system particularly at the edge of the 
cloud, where the TF solution slightly overestimates the 
size of the outer fully paired region. For smaller sys- 
tems we found that the agreement deteriorates. Indeed, 
at small N a partially polarized state extends until the 
edge of the cloud even at intermediate coupling and small 
spin imbalance, in contrast with the TF approximation, 
which predicts a phase separation with fully paired wings. 

As pointed out in a number of previous work a 28 i 29 ' 37 i 38 
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FIG. 10: (Color online) Local polarization profile for N = 200, 
P — 4%, g = 20, at various temperatures T measured in 
units of Tujj z . For comparison, the central density in the trap 
corresponds to an energy cf = h 2 /2rnk F ~ 272. The inset 
magnifies the region of the cloud around the polarization dip, 
where the pairing is the strongest. At this stronger coupling, 
the fully paired region is more stable, and starts breaking at 
a higher temperature (above T = 3.5 is this case). The filled 
triangles indicate the polarization at the cloud radius which 
contains 99.9% of atoms. 

the pair momentum distribution reveals the signature of 
the FFLO inner shell with two symmetric peaks centered 
at — kp\. It could be measured in the experiment 

by quickly removing the confining lattice and sweeping 
the magnetic field to the BEC side of the Feshbach reso- 
nance. This would bind the pairs into molecules, whose 
momentum distribution could be measured by means of 
time-of-flight imaging techniques^ The tails of the one- 
body momentum distribution strongly depend on the in- 
teraction strength and decay as 1/fc 4 , due to the non ana- 
lyticity of the exact wave function as a consequence of the 
delta function interaction. This could be another inter- 
esting quantity to measure in the experiment, although 
it does not carry information on the actual pairing in the 
system, but only on the effective strength. 

The finite temperature PIMC calculations done with 
N = 200, P = 4%, and three different couplings (g = 8, 
16, and 20) show a crossover from a phase separated state 
to a partially polarized state at T ~ 0.065 T F for g = 8, 
and T ~ 0.10 T F for g = 20. For the strongest coupling- 



analyzed with this method (g = 20), the fully paired re- 
gion disappears already at T ~ 0.035 T F , leaving a dip in 
the polarization density profile which becomes shallower 
and shallower as the temperature increases. This quite 
low temperature (the experiments will access tempera- 
tures down to T ~ 0.05 T F ) puts some limitation on the 
visibility of the fully paired phase at weak-intermediate 
coupling. One has to go to stronger coupling to make the 
pairs in the outer shell stable up to higher temperatures, 
being the binding energy proportional to g 2 . 

We discussed the ergodicity issues in the Monte Carlo 
sampling at strong coupling, and we identified the prob- 
lem with the difficulty of unpaired fermions to pass 
through strongly bound pairs. The same slowing down 
in the equilibration might occur also in the experiment 
for ultra thin traps (Nhuj z <C Tiu>± and ksT <C Twj±_). 
On the other hand, at weak coupling there are few par- 
ticles involved in the fully paired region, which is very 
narrow, and the signal coming from this phase should be 
very difficult to detect in the experiment. We suggest 
that a good spot to detect the transition from the fully 
paired to the fully polarized state is for effective couplings 
g/VN in the range 3 — 5, where the system features a 
fully paired region up to 10% of spin imbalance, the size 
of the two phase separated states is quite large, and the 
equilibration is not hard to reach. 

The information coming from this paper could be com- 
bined with the recent work by Parish et aL,— who studied 
the quasi one dimensional effects coming from the array 
of tubes in the transverse direction by means of a tight- 
binding model, in order to find out the parameters of the 
Hamiltonian which maximize the visibility of the FFLO 
state. In principle, one can extend the present work by 
explicitly including an array of tubes in the QMC simu- 
lations to study quasi one dimensional effects present in 
the experimental setup. 
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